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Abstract 

In this paper we consider the special case where a signal x G is known to vanish 
outside a support interval of length m < N. If the support length m of x or a good 
bound of it is a-priori known we derive a sublinear deterministic algorithm to compute 
X from its discrete Fourier transform x G C'^. In case of exact Fourier measurements 
we require only O(mlogm) arithmetical operations. For noisy measurements, we 
propose a stable 0{mlogN) algorithm. 
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1 Introduction 

It is well-known that FFT algorithms for the computation of the discrete Fourier 
transform of length N require 0{N log N) arithmetical operations. However, if the 
resulting vector is a-priori known to be sparse, i.e., contains only a small number of 
non-zero components, the question arises, whether we can do this computation in an 
even faster way. In this paper we derive a new deterministic sparse FFT algorithm 
for signals x G that are a-priori known to vanish outside a support interval of 
length m < N. 

Related work. In recent years many sublinear algorithms for the sparse Fourier 
transform have been proposed that focus on (approximately) m-sparse vectors or 
vectors being norm bounded. 

Most of the these methods are randomized poly-logarithmic sparse Fourier algo¬ 
rithms achieving e.g. a complexity of 0{mlogN) [5] for m-sparse signals, or even 
0{mlogm), see e.g. [HIUS]. However, these algorithms succeed to compute the cor¬ 
rect result only with constant probability being smaller than 1. Moreover, there is 
no efficient method to check the correctness of the result. Another drawback is that 
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samples need to be drawn randomly, this is significantly more complicated to realize 
by hardware. Regarding the numerical results, runtimes start to be more efficient 
than standard FFT for m = 50 and N > 2^”^, see mm for most of the considered 
algorithms, and for m = 50 and > 2^^ for the algorithm in [6]. Another exper¬ 
iment shows for fixed N = 2^^ that several sparse FFT algorithms start to pay off 
for m < 2200 while the algorithm by [U] is more efficient for m < 2 ^^, mi. Newest 
algorithms are even faster, but the probability to provide the exact result decreases 
with larger m. For a nice overview of the techniques of randomized sparse Fourier 
transforms we refer to [3]. 

Deterministic sparse FFT algorithms proposed e.g. in p n [g ini El] are appli¬ 
cable to (approximately) sparse or so-called norm bounded signals. Iwen [9l Enj con¬ 
sidered sparse m-term approximations of the discrete Fourier transform and achieved 
an i^jy/m error bound with 0(m^ log^ N) operations. Similar runtime is achieved 
for the deterministic algorithm in where accessibility not only to the usual signal 
samples but also to time-shifted samples of the input function is assumed. The deter¬ 
ministic algorithm in [ 2 | for signals being norm bounded by \/m, the evaluation of 6- 
approximations of all r-signihcant frequencies, has polynomial costs in log A^, m, Ijr 
and 1/5. 

Finally, we mention sparse FFT algorithms based on Prony’s method that are 
completely deterministic and can be operated with 0{m^) operations (independent 
of the signal size N), see e.g. [TIES], but usually with an unstable numerical behav¬ 
ior. A better numerical performance can be achieved by randomization leading to a 
complexity m^/^log^ A^, see [Tj. 

Compared to the approaches above, we restrict ourselves to signal vectors possess¬ 
ing a small support interval. Vectors with small support appear in different applica¬ 
tions as e.g. in X-ray microscopy, where compact support is a frequently used a-priori 
condition in phase retrieval, as well as in computer tomography reconstructions. 

Notations. First, we hx some notations. For a vector x G of length N , 
let its discrete Fourier transform be defined by x = Ftvx, where the Fourier matrix 
Ftv G is given by 

Fat := , ajN:=e~. 

\ ) j,k=0 

Let the support length m = |suppx| of x G be defined as the minimal integer 
m for which there exists a. p, G {0,...,Ai — 1} such that the components of x 
vanish for al\ k ^ I := {{p + r)modN, r = 0,..., m — 1}. The index set I is called 
support interval of x. Obviously, the support length m of x is an upper bound for 
its sparsity ||x||o, i.e., the number of nonzero components of x, since there may be 
zero components of x inside the support interval I. However, we always have 7 ^ 0 
and / 0. Observe that the support length m of a vector x G is always 

uniquely defined while the support interval itself resp. the first support index p needs 
not to be unique. For example, considering the vector x = {xk)^SQ G of even 
length N given by xi = XAr/ 2 - 1-1 = 1 and = 0 for A: G {0,..., V — 1} \ {1, N/2 + 1}, 
we find for x the support length m = V/2-1-1 while the support interval can be chosen 
either as {1,..., N/2 + 1} or as {N/2 -|- 1,..., V — 1,0,1}. However, if m < y, then 
the support interval and hence also the hrst support index p are uniquely determined. 
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Our contribution. In this paper, we will derive new deterministic algorithms 
to reconstruct a vector x G with support length m < N from its discrete Fourier 
transform x G C'^. In case of exact data, we will use less than 4m Fourier samples 
and 0(mlog2m) arithmetical operations to recover x. Thus, the algorithm already 
starts to pay off for m < 

In case of noisy Fourier measurements, we will use 0(mlog2 A^) arithmetical op¬ 
erations. More exactly, we will employ one or several FFTs of size less than 4m as 
well as the computation of [log2 — 1 scalar products of length m to recover the 
full vector in a stable way. 

In both cases, for exact as well as for noisy measurements, the algorithms are 
based on the idea that for N = 2"^ the nonzero components of x can already be 
computed by evaluating a periodization of x with length 2^ > m. Thus, for the 
complete reconstruction of x, we only need to compute in a second step the correct 
support interval resp. the first support index of x. 

Organization of the paper. In Section 2, we recall some properties of the 
discrete Fourier transform that will be used in our approach. Section 3 is devoted to 
the sparse FFT algorithm for m-sparse vectors in case of exact Fourier measurements. 
Section 4 considers a more stable variant of the different steps of the algorithm that 
make the method robust in presence of noise, and even improves the accuracy of the 
resulting vector in comparison to the usual FFT method while using only 0(mlog A^) 
arithmetical operations. Finally, we present our numerical experiments in Section 5. 


2 Preliminaries 

Throughout the paper let N := 2^ with some J > 0. 
We consider the periodizations G of x, 


Mi) _ _ 

~ A=o ~ 



2^-1 


k=0 


( 2 . 1 ) 


for j = 0,..., J. Obviously, x^*^) = Ylk=o is the sum of all components of x, x^^i = 
^ and = x. We recall the following relationship for 

the discrete Fourier transform of the vectors x^^\ j = 0,..., J, in terms of x. 


Lemma 2.1 For the vectors x^-^i G j = 0,J, in (2.1), we have the discrete 
Fourier transform 

-Fn xi-^i — (r T 

X .— ^23^ — \^2J-^k)k=Q ’ 


where x = = F^rx is the Fourier transform of x £ . 

Proof: By definition, we find for the components of x^-^i 


2^-1 


2J-1 


^ =Y. Y. ^r+2l,‘ 


xr' : = 


IUIn 


r=0 

N-l 


r =0 £=0 


2'' ^nk 

/ ^ Xn — •^2'’-3ki 


A: = 0, 


, 2 ^- 1 . 


n =0 
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Thus the assertion follows. ■ 


Lemma 2.1 tells us that for a given vector x the Fourier transform of the peri- 
odized vectors needs not to be computed but is just obtained by picking suitable 
components of x. Further, we want to make use of the following simple observation. 


Lemma 2.2 Let x = G C^, N = 2-^, and let y = {yk)k=o ^ ^e its 

shifted version with components yk ■= a:(fc+ 2 iy)modAr; k = 0,...,A^ — 1, for some 
j G {0,..., J — 1} and v G {0,..., 2^~^ — 1}. Then the components of the Fourier 
transformed vectors x = F^rx, y = F j^y satisfy 

In particular, for j = J — I and n = 1 we have yk = a^(fc+Af/ 2 )modiV; A: = 0,..., — 1, 

with 

^ ^ ^ ^ , N 

X2k = ?/2fc, X2k+1 = —y2fc+l, A: = 0, . . . , ^-1. 

Proof: With x = Fatx = {xi)^S^ and y = FAry = we obtain 

7V-1 N-l N-l 

Ik sr^ Ik sr^ l{k-2iu) -h, ^ 

yi= 2^ yk^N = X{k+2Jy)raodN = 2^XkUjf^ ’ = Xi- 

k=0 k=0 k=0 

For j = J — 1 and n = 1 the assertion follows with = uj~^ = (—1)^ ■ 


3 Reconstruction of x from exact Fourier data 

We assume that the Fourier transformed vector x = F^rx is given, and that the 
support length m of x or an upper bound m of it is known a-priori. Now, we apply 
the following simple procedure to recover x. 

Let < m < 2^. Then, by Lemma 1, = ix2J-{L+i)2)k=o 

Fourier transform of In a first step, we compute using inverse FFT of 

length 2^~^^. 

Since |suppx| = m < 2^ by assumption, it follows that has already the 

same support length as x since for each k G {0,..., 2^~^^ — 1} the sum in 


r\J — L — l _ 1 

(L+i) _ sr-^ 

Xk ~ / ^ ®fc-|-2^+l£ 

1=0 


(3.1) 


contains at most one nonvanishing term. Moreover, also the support itself, or equiv¬ 
alently the first support index n = of x(^+^), is uniquely determined. 

Thus, in order to recover x from we only need to compute the correct first 

support index of x, such that the components of x are determined by 


^ +k)raod N 


From (3.1) we know that 
{0 ,l,...,2^-^-i-l}. 


^(/x(^+l)+A:)mod2^+i 

0 

is of the form 


A: = 0,..., m — 1, 
k = m,... ,N — 1. 

= -|- 2^~^^v for some v G 
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Theorem 3.1 Let x G , N = 2'^, have support length m (or a support length 
bounded by m) with < m < 2^. For L < J — 1, let be the 2^+^- 

periodization of x. Then x can be uniquely recovered from and one nonzero 

component of the vector {x2k+i)^=Q 

Proof: Let have the support interval {(/x(^'’'^)+r)mod2'^'’“^, r = 0,..., m—1}. 

Since m < 2^, this support interval is uniquely determined. Further, we know that 
by construction the hrst index of the support interval of x has the form = 
^(^+1) _|_2^+izy for some G {0,..., 2'^“'^“^ — 1}. We now consider the vector G 
that is given by the components 


7/0 

«^(i+i)+fc 


(L+l) 

^(/7(^+i) +fc)mod 2 ^+i 

0 


A: = 0,..., m — 1, 
k = m,..., N — 1. 


The vector is obtained as one possible vector in with support length m possess¬ 
ing the 2 ^^^-periodization Obviously the first index of the support interval 

of uO is Therefore, all further vectors in with a support length m and 

possessing the 2 ^+^-periodization can be written as a shifted version of of 

the form 


:= iuk)k=d with 


~ '“(fe-|- 2 i+li/)modAf’ k — 0, . . . ,N 1, 


for some z^g{ 1,...,2'^ ^ ^ — 1}. Thus, the desired vector x is contained in the set 
: u = 0,...,2'^-^-^ - 1 }. 


It remains to find u such that x = u^. From Lemma 2.2 it follows 
nents of the Fourier transform that uf = u^*. 

Choosing now an odd-indexed nonzero Fourier value X 2 fco+i the 
we can compare it to the corresponding component of 


for the compo- 
given vector x, 


m—1 

^0 (-t'+l) (//('^+l)+r)(2A:o+l) 

^2fco+l = 2^ ^(^(^+l)+r)mod2^+l 
r=0 


and obtain the correct value v from = 5l2fco+i/^2fco+i- remark that 

zz G {0,..., — 1} is indeed uniquely determined by since gcd(2fco -b 

1,2'^“^“^) = 1, i.e., 2kQ -\-1 and 2'^~^~^ are mutually prime. Finally, the vector x is 
obtained as x = u^. ■ 

We summarize the algorithm to evaluate x from x for exact Fourier data as follows. 


Algorithm 3.2 (Sparse FFT for veetors with small support for exact Fourier data) 
Input: X G C'^, N = 2“^, |suppx| < m < N. 

• Compute L sueh that 2^~^ < m < 2^, i.e., L := |'log 2 m]. 

• If L = J or L = J — 1 , compute x = F^^x using an FFT algorithm of length 

N. 

• IfL<J-l: 
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1. Choose a := (5?2-^-(i+i)fe)fc=o ^ compute 

using an FFT algorithm for the inverse discrete Fourier transform of length 
2 ^+ 1 ^ 

2. Determine the first support index G {0,..., — 1} of 

such that / 0 and = 0 for k ^ + r)mod 2 -^+^, r = 

0,..., m — 1}. 

3. Choose a Fourier component X2ko+i / 0 o/ x and compute the sum 


'O 

U2ko+i 


SL+l) 


m—1 

®(At(i+l)+£)mod2i+l 

e=o 


{2ko+l){ti‘-^+^''+i) 


N 


4- Compute the quotient b := X 2 fco+i/^fco+i lhat is by construction of the form 
b = for somep G {0,..., — , and find G {0,..., 2'^“^“^ — 

1 } such that {2ko + l)v = pm.od2^~^~^. 

5. Set := + 2 '^+^z/, and x := {xk)^Fn with entries 


^(/^(■^)+£)modAf ■ 


X 


(L+l) 

(^(L+l)_l-£)mod 2-^+1 


0 


= 0,..., m — 1, 
= m,..., N — 1. 


Output: X. 


We simply observe that the proposed algorithm has an arithmetical complexity 
of Cl(mlogm). In step 1 we employ an FFT algorithm of length 2^^^ < 4m having 
this complexity. All other steps can be performed with 0{m) or less operations. 
Particularly, the support interval (resp. the first support index of in 

step 2 can e.g. be found by computing the local energies := Y^^=k~^ I®£^d 2 ^+il^ 
for A; = 0,..., 2 '^+^ — 1 , and taking ;= argmax;;, Ck- Here, cq can be computed 

by at most 0{m) operations, and all further energies by using the recursion Ck+i = 
efc ~ \xk\‘^ + \xk+m\^^ where we need to keep in mind that there are only m nonzero 
entries in xH+i), 

The complete algorithm requires less then 4m Fourier values for the overall eval¬ 
uation. 


Let us give some further remarks on Algorithm 3.2 


Remark 3.3 

It is always possible to choose a nonzero Fourier component of the vector 
{x2k+i)kLo This can be seen as follows. Let suppjt. = {p^T ^ ^Op-iin.odN,p^T-\- 
m — ImodA"} be the support of x, where m < 2^ < N/A. Then the trigonometric 
polynomial 


p{u) = 


N-l 

E 

k=0 


Xk e 


— \u)k 




(J) 


m—1 

E 

£=0 


^+£) mod W 


iuji 


(3.2) 


is real, non-negative and of degree < m — 1. Hence it possesses at most m — 1 
pairwise different zeros, where all zeros are at least double zeros. Observing now that 
= P 'xe can conclude that not all X 2 k+i, k = d,... ,N/2 — 1, can be zeros 

of p{u)) since N/2 > 2m. 
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Remark 3.4 

As shown in the remark above, there can be at most m — 1 zero components in the 
vector {x2k+i)kL^ However, for a stable computat ion of the correct first support 
index the Fourier value X 2 fco+i used in Algorithm 3.2 should have large modulus. 
This may be ensured by taking 

/co := argmax{|x 2 fe+i|^ : k = 0,.. .,N/2 - 1}. 

Unfortunately, this procedure is quite costly. Instead, we may compute 

:= argmax{|x2J-(i+i)fc|^ : k = 0,..., 2^+^ - 1}. 


Then we have for the trigonometric polynomial in (3.2), 


\x 




2J-(i+ 


and it is likely that p I ^l+i — ) is close to the global maximum ||p||oo of p{uj). In 

step 3 of Algorithm |3.2| we may now choose the one neighboring Fourier value with 
maximal modulus, i.e., ezi/ier 

Ifp{uj) attains its global maximum at some point uiq G 


then it follows that 


UJo 


2Tik, 


{L+l) 


27r(4^+^^-l/2) 27r(4^+^^ + l/2) 


2L+1 


2L+1 


2 L + I 


< above choice o/x2fco+i with ko = 


2 “^ ^ or fco = 2 '^ ^ — 1, we can further assume that 

2 iffT, und applying the Lemma of Steckin (see e.g. we find 


Wo 


27r{2/i:o+l) 

N 


< 


\X2ko+l\ =P 


2 _ ^ / 27 r( 2 fco+l)\ 


V 


N 




, T(m—1)\ ^ r\ 

, cos ( -^r+r^ I > 0. 


4 Reconstruction of x from noisy Fourier data 

Let us now assume that the given Fourier data are noisy, i.e., they are perturbed by 
uniform noise, 

yk = Xk + ek (4.1) 

with |efc| < 6. Similarly as before, we want to reconstruct x from the noisy Fourier 
vector y using the a-priori knowledge that x has a support interval of length m < N. 


For that purpose, we propose a stabilized variant of Algorithm 3.2 The stabiliza¬ 
tion regards the following essential steps in the algorithm, namely 

( 1 ) the correct determination of the support interval of x(^+^\ 

( 2 ) the correct determination of the support interval resp. the first support index 
of the desired vector x, and 

(3) the evaluation of the nonzero components of x within the support interval that 
may be improved by employing more Fourier values y than in the case of exact data. 


(1) Stable determination of the snpport interval of 

For that purpose we use the following observation. 
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Theorem 4.1 Let x G , N = 2'^, have support length m (or a support length 
bounded by m) with < m < 2^. For L < J — 1, let 

be the 2^~^^-periodization of x and x(^+^) = (J 2 ''-i-ifc)fc=o Fourier transform. 

Then, for each shifted vector of the form 

iW := (x2J-i-ifc+K)Lo ’ K = 0,..., - 1, 

the inverse Fourier transform satisfies 

Proof: By definition, we obtain for the components of 

2 i+l_i 


(k) 

Zi = 


1 


2^+1 

1 

2-^+1 


X2-^-C-ik+K, ^2^+1 

fc =0 

2 ^+1-! AT-l 

r{2-’-Z-^k+K) ^ -ki 
2 L + I 


Y1 


OJ. 


1 


k=0 r=0 

N-l 


2L+i_i 


2 L+I 


E E 


UJ. 


k{r—i) 

2 L +1 


r=0 
2J-Z-i_i 

j=0 


k=0 


(£+ 2 i + lj)K _ 
^AT 


1-1 


CJ'at ■" — Wjv ^ ^ 3l£_|_2i+ljW 

j=0 


2 ^ + ljK 
IV 


Since |suppx| = m < for each £ the above sum contains only one value, thus 


4^)1 


(i+i)i 


\zl '\ = |x^ |. - 

Obviously, we have by definition. We observe that all vectors 

are constructed from different Fourier components of x. This circumstance will be 
used to stabilize the algorithm. Applying the noisy measureme nts ijk in (4.1) instead 

of Xk, let 'z := {y 2 J-v-ik+K)'kJo before in Algorithm 3.2, we compute in a 

first step the vector = y(^+i) from the noisy measurements {y2J-(r+i-)k)‘k=o 
In order to determine the support interval of from the vector z^*^^ we consider 

the energies 


2 l( 0 ) .- 
^k ■ — 


m-\-k—l 

\y. 

£=k 


(L+l) 
£mod2^+l I 


m+k—1 

|2 _ lAO) 


E 

i=k 


Z 


t mod 2^+11 ’ 


A: = 0,..., 2^+1 - 1 


and take ■— argmaxe^^^ as the first estimate for For higher noise 


levels, we stabilize the computation of as follows: We compute also the vector 

J-L-2\ 

) ^ 2 ^ + l-l 


J <2 _ ^^2 

z(^ ) by applying the inverse FFT to z 


the energies 


m+/c —1 

t mod2^+11 ’ 
'=k 


^'k' ■= Y1 


= {(j2J-z-^2k+i))k=o ’ compute 

A: = 0, ...,2^+^ -1, 







then this index is taken 


and take '■= argmax |+e^^'’) 


1I 




as the first support index of Otherwise, we compute a further vector, e.g., 

2J-L-3(^4k^i'j)‘lJQ evaluate the corresponding energies and 






:= argmax |(e 


_l_ _|_ 

sV^fc 


etc. 


(2) Stable determination of the support interval of x. 

In order to stabilize the computation of the first support index of the complete 
vector X = x^*^) resp. the shift v such that we employ an 

iterative procedure, where at each iteration step, the support interval of the periodized 
vector x(-^"’“^\ j = L + 1,... J — 1, is computed from the support interval of the vector 
x^-?) of half length. 

At each iteration step we apply the following procedure. Observing that x^-^'*'^) 
has the same support length m as and using the relation 


X 


(i+i) I Jj+^) 


+ X 


k+2i 


= X 


(j) 


k = 0,, 


, 2 ^- 1 , 


it is sufficient to check whether or + 2-^, i.e., whether the 

first support index of x^-^^ is equal to the first support index of or 

if the support has to be shifted by 2^. We illustrate the two cases in detail in Figure 

m 


First case: 

xO') ___ 


X 


(i+1) 


X 


(j) 


X 


O'+i) 


Second case: +2K 

x(i) ___ _ 


X 


(i+i) 


X 


(f) 


X 


(i+1) 


Figure 1: 

Possible support change in one iteration step. 


Generally, in order to recover x^-^'*'^) from we apply the following procedure. 

Theorem 4.2 Let x G C^, N = 2'^, have support of length m (or a support length 
hounded by m) with 2^~^ < m < 2^. Fu rther, let L + l<j<J — 1, be the 
periodizations of x = x^'^^ as given in (2.1). Then, for each j = L + 1,..., J — 1, the 


vector can be uniquely recovered from x^^'i and one nonzero component of the 

vector of Fourier values y ;= ix 2 j-u+i)( 2 k+i))‘k=o ■ 

Proof: Let x^^'l be the given vector of length 2^ with support suppx^-^^ = + 

r)mod2'^, r = 0,..., n — 1}. In order to determine we only need to decide 
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whether or + 2^, or equivalently, whether is given 

by 


™(i+l) _ f ®(/iO)+£)mod 2 J ^ — 0 ,...,n 1, 

I 0 else, 

or by 

(i+l) _ f ®(^(j)+€)mod 2 J ^ = 0, ...,n — 1, 

(/xO)+ 2 J+(?) mod 2 J +1 I Q else. 


(4.2) 


(4.3) 


The two possible solutions only differ by a shift of all components by 2T For 

sim plicity let us denote the two solutions by and . Then, according to Lemma 
the Fourier transformed vectors = (m®)^'’Jq~^ and satisfy 


2.2 


the relationship 


40) 


41 ) 


^2fc+l ~ ^2fc+l’ 


A: = 0,..., 2^' - 1. 


Using now one nonzero Fourier value i?2-^-(j+i)(2fco+i) ~ ^ 2 ko+i^ determine 

the correct vector x^^^^. For that purpose, we just compute 


u. 


ra—1 

( 0 ) _ Jj) 


2fco+l 


E 

1=0 


®(/i(j)+£)mod2j‘^2J+l 


(2fco+l)4+MO)) 


and compare it to the Fourier value )?2"'-0+i)(2fco+i)- l^2feo+i ~ ^2-^-o+i)(2fco+i)l 
l^2fco+i + ^2‘^-(j+i)(2fco+i)l’ then x^+i) = and = ^U)^ otherwise, we find 

x4+i) = and ;ub'+i) = ^(i) + 24 ■ 

(3) Evaluation of the nonzero components of x. 

Finally, we can use the vectors G ^ computed in step (1) also for a more 
exact evaluation of the nonzero components Xk of x as follows. First, we employ our 


investigations in Theorem 4.1 and observe that 

^ ^ ^{L+l) ^ 

Using the equation we can reformulate this as 


+('^) _ , (rr , + 

^£mod2L+i ~1^4+2'^+1+modAf h 


~{k) 

Z, 


+k) 


k = 0,... ,m — 1. 


‘'+(r+l)+fc) mod 2^+1 mod V 

Therefore, if we have to compute the support entries a:+(J)+fc) mod V) k = 0,... ,m— 1, 

) N-1 
k=0 


from the noisy measurements (Sk)k=o take the average 


^+k) mod V 


Q-i- \ ^ ^ (/4,(-^+i)+/ii) mod 2^+1 ’ 

r=0 


A: = 0 ,..., m — 1 , 


where i? + 1 is the number of the vectors z('") = 'P^l+i{y2J-L-ik+K)k=Q that we 
want to involve. 

The complete algorithm to compute x from its Fourier transform by a fast sparse 
FFT algorithm can now be summarized as follows. 
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Algorithm 4.3 (Sparse FFT for vectors with small support for noisy Fourier data) 
Input: noisy measurement vector y G , N = 2'^, |suppx| < m < N. 

• Compute L such that < m < 2^, i.e., L := |'log 2 m]. 

• If L = J or L = J — 1, compute x = by inverse FFT. 

• IfL<J-l: 

1. Choose 'z ^ := y(^+^) = (y 2 -^-(i+i)fc)fc =0 o,nd compute = y('^+^) := 

using an FFT algorithm for the inverse discrete Fourier trans¬ 
form. 

2. Determine the first support index G {0,..., 2^~^^ — 1} using 

the following iteration: 

• Compute the energies 


^ 0 ) — 


E 

e=k 


-^ 0 ) ,2 

£ mod 2^+11 ’ 


Z, 


and compute '■= argmaxe^°\ 


k = 0, - 1 


Compute ^ by IFFT from (y2''-i-2(2fc+i))fc=o determine 


m+k—1 

Z_^ I mod 2^+1 1 
e=k 


k = 0, - 1 


and take '■= argmax | 

k 

• Set j := 0. 

WhUe 

proceed by computing for a further k G {1,..., 2'^~^~^ “ 1} \ {2'^“^“^} 
the vector by IFFT from (y 2 "'-^-ifc+K))fc=o energies 


m+fc—1 

2(i+2) |~(«) ,2 

f mod 2^+11 ’ 
=k 


e'k ' ' ■= Y. 


k = 0,..., 2^+^ - 1 


and take /ij+ 2 ^^ •= 

k 

Set := and j := j + 1. 

End (while). 

S'et x(^+^) = 


X 


(i + l) 

(fc+/^(-^+l)) mod 2-^+1 


-J mod 2^+1 ^ = 

1 0 k = m,..., 2^+^ — 1. 


3. For j = L + 1,..., J — 1 

Choose a Fourier component y 2 J-{j+^)( 2 ko+i) ~ ^ ^ compute 


m—1 


®i+i •” 'Yj 

£=0 


mod 2-^+1 


^2J+l 
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If l«i+i - + y^C+il’ x(^+^) := 

i^^k^^"’)‘k=o~^ wit/i entries 

^(i+i) ^ f ^(;^ 0 )+£)mod 2 i ^ = 0, ...,n-l, 

\ 0 else. 


If\aj+i-y^^i^^^\ ^ l“i+i + yEt+il’ +2-’' andx('^'+^) := 

with entries 


.j. 0 '+i) 

(ajO) +2J +£)mod 2J+1 


^(/xO)+£)mod 2 J 
0 else. 


4- Assuming that G , r = 0,. 

step 2, we compute 


, B, have been evaluated already in 


B 


i = fj. 

Output: = X 




^(£+2^+1 1 /) mod N 
^(■f'+l) _|_ jTT, _ x_ 


^ A^r) , -l^rie+^^ + '^u) 

r=0 


Let us shortly summarize the arithmetical complexity of the algorithm. Step 1 
and 2 together require two or more inverse FFTs of length < An, i.e. O(mlogm) 
arithmetical operations. The number of involved vectors depends on the noise 
level. The evaluation of energies can be done at each iteration step by 0(m) opera¬ 
tions. 

In step 3, J — (L + 1) = log 2 N — |'log 2 m] — 1 < log 2 (A^/ni) iterations are needed, 
where at each iteration a scalar product of length m has to be computed and to 
be compared to a given Fourier value. To find a nonzero Fourier value needs less 


than m comparisons, see also Remark 3.4 Hence only 0(m log m + mlog(N/m)) = 


0{m log 2 N) arithmetical operations are needed to recover x in the case of noisy data. 


Remark 4.4 

For the computation of x it is sufficient to compute as well as the first 

support index pf'^\ where j = L + 1,J — 1, is iteratively computed from 

pfA. The values Oj+i can be obtained directly from x(^+^) and . Hence , th ere is 
no reason to compute the intermediate vectors x^-^'*'^) in step 3 of Algorithm 4-3, they 
are only given for better illustration. 


5 Numerical results 


In this section, we discuss the behavior of the algorithm for noisy input data. We 
reconstruct randomly generated vectors x from disturbed Fourier data y = x -|- £, 
where we assume uniform noise s = {ek)^TQ with \ek\ < 5 at different noise levels. 
As a noise measure, we use the SNR value 


SNR = 20 • log 


10 


|X||2 

l^lb' 
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Let us first illustrate the algorithm for a vector x of length N = 2^ = 256 and support 
length m = 6, i.e., we have J = 8 and L = 3. The nonzero entries of x are xio 5 = 8, 
3^107 = xio8 = —5 and xno = 2. We add a noise vector e of SNR = 20 to 
X and reconstruct x from y = x + £. For the noise vector e in this example, it 
holds that ||£||oo = 1-721 and ||£||i/256 = 0.952. In Figure]^ we present x, and we 
compare the reconstruction results of our algorithm to the result of the inverse Fourier 
transform applied to y. The reconstruction x' computed by our deterministic sparse 
FFT algorithm has nonzero components = 7.884+0.131 i, x'^Qg = —0.070—0.149 i, 


no? 


= -3.100 + 0.171 i, x'los = -4.955 + 0.094 i, = -0.105 + 0.022 i 




1.879 — 0.076 i, and an error ||x —x'||2/256 = 0.00146 whereas the error by the inverse 
FFT is ||x — F25gy||2/256 = 0.00395. In this example, the reconstruction by our 
algorithm requires no additional vectors for the determination of in step 2 of 


Algorithm 4.3, i.e., we only use the two vectors and l 2 ) _ ^(s) q£ 


2^~^^ = 16 in this step. Thus, we have taken 32 + 4 = 36 Fourier values to recover x. 
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Figure 2: 

(a) Original vector x of length N = 256; (b) Reconstruction of x using the sparse 
FFT Algorithm 4.3 (c) Reconstruction of x using the inverse FFT. 


In Figure]^ we show the reconstruction error using the sparse FFT Algorithm |4.3| 
for vectors x of length N = 2^^ and support length m = 50 where the vectors are 
randomly generated with |Re(xfc)| < 10, |Im(xfc)| < 10 for k in the support interval. 
We consider noisy Fourier data with SNR values between 0 and 50 and compute 
the reconstruction x! for 100 vectors x at different noise levels. The quality of the 
reconstructed vectors x! is evaluated by computing the norm ||x — x'|| 2 /A^. Figure 
shows the average error norm over all 100 considered vectors. We compare the 
reconstruction results of our algorithm to the results of an inverse Fourier transform 
applied to the noisy vectors y. For noise levels SNR < 10, our algorithm made use of 
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additional vectors in step 2 of Algorithm |4.3| in order to improve the identification 
of the hrst support index of Here, the algorithm evaluated at most five 

additional vectors (for SNR = 0), hence a total number of up to seven vectors 
has been used to determine in some cases. 



Average reconstruction error ||x — x'|| 2 /iV for different levels of uniform noise, 
comparing our deterministic sparse FFT algorithm and inverse FFT. 


Determining the first index of the support of x (i.e., finding the correct /r = 
is one of the crucial points of the algorithm for noisy input data. If is not 

identified correctly, the correct support interval cannot be found anymore, even if all 
shifts in step 3 of Algorithm 4.3 are correct. 

In a third experiment we again take randomly generated vectors x of length N = 
2^^ with support length m = 50 or m = 2^®, where |Re(xfc)| < 10, |Im(xfc)| < 10 for 
k in the support interval. 


SNR 

N = 2^^ m = 50 

N = 2^^ m = 2^® 

correctly 
identified fi 

ll^lloo 

M\i/N 

correctly 
identified fi 

II ^ II oo 

M\i/N 

0 

86% 

68.367 

37.002 

78% 

4927.294 

2666.891 

5 

97% 

37.799 

20.458 

93% 

2911.275 

1575.695 

10 

99% 

22.262 

12.049 

97% 

1365.686 

739.186 

15 

100% 

12.559 

6.798 

100% 

841.737 

455.585 

20 

100% 

6.751 

3.654 

100% 

483.223 

261.542 

25 

100% 

3.796 

2.055 

100% 

261.897 

141.752 

30 

100% 

2.147 

1.162 

100% 

144.593 

78.259 

35 

100% 

1.242 

0.672 

100% 

84.374 

45.665 

40 

100% 

0.668 

0.362 

100% 

47.666 

25.800 


Table 1: 


Percentage of correctly identified /r and average norm of noise in 100 randomly chosen 
vectors for different noise levels, dependent on length N and support length m of x. 
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Table shows the number of cases in which the first support index fj, has been 
determined correctly by Algorithm |4.3| for 100 randomly chosen vectors for each noise 
level. Additionally, we give an average for the norm of the noise vectors s at each 
noise level. 

The results show that the algorithm succeeds for very short support intervals as 
well as for long support intervals compared to the full vector length. In cases where 
/i could not be determined correctly, the error originates from step 2 of the algorithm 
where has to be determined. However, the deviation from the correct support 

was small (< 6) in any case. 

We can conclude that even for high noise level, the support of the reconstructed 
vector X is correctly found in most of the cases. The support is always correct for 
noise levels with SNR > 15. Table also shows that the absolute noise at each 
component can be considerably large. It is essentially larger for vectors with larger 
support since in this case also the signal energy grows accordingly. 
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